Changes to the gut microbiota of a wild juvenile passerine in a multidimensional urban mosaic

Urbanisation is a major anthropogenic perturbation presenting novel ecological and evolutionary challenges to wild populations. Symbiotic microorganisms residing in the gastrointestinal tracts (gut) of vertebrates have mutual connections with host physiology and respond quickly to environmental alterations. However, the impact of anthropogenic changes and urbanisation on the gut microbiota remains poorly understood, especially in early development. To address this knowledge gap, we characterised the gut microbiota of juvenile great tits (Parus major) reared in artificial nestboxes and in natural cavities in an urban mosaic, employing two distinct frameworks characterising the urban space. Microbial diversity was influenced by cavity type. Alpha diversity was affected by the amount of impervious surface surrounding the breeding location, and positively correlated with tree cover density. Community composition differed between urban and rural sites: these alterations covaried with sound pollution and distance to the city centre. Overall, the microbial communities reflect and are possibly influenced by the heterogeneous environmental modifications that are typical of the urban space. Strikingly, the choice of framework and environmental variables characterising the urban space can influence the outcomes of such ecological studies. Our results open new perspectives to investigate the impact of microbial symbionts on the adaptive capacity of their hosts.


Results
We collected faecal samples from 15-day old great tit nestlings hatched in nine different sites, located within and outside of the capital city of Warsaw in Poland (Fig. 1). We characterised gut microbial communities using 16s ribosomal RNA (rRNA) gene sequencing. After quality-filtering, the resulting dataset consisted of 82 Table 1). Similarly, microbial composition did not differ between nestlings originating from nestboxes and natural cavities (PERMANOVA: Bray-Curtis: pseudo-F = 0.25, p = 1; weighted UniFrac: pseudo-F = 0.53, p = 0.92).
The impact of urbanization on microbial community diversity and composition. To understand the impact of urbanization on the gut microbiota, we defined the urban space using two different frameworks. Under the first framework of urbanisation, each nest was categorised depending on whether it was located within the administrative limits of the city (urban) or outside of these limits (rural, Fig. 1). This aligns our findings with frequent partitioning of phenotypic 20 and genetic variation 21 . However, the actual surface use patterns measured by the percentage of Impervious Surface Area around each nest location (ISA) show substantial heterogeneity within urban and rural sites (Fig. 3a, Supplementary Table 2). Furthermore, ISA correlated with all the environmental parameters used in the study (Fig. 3b, Supplementary Table 3). Taken together, ISA can be used as an urban metric to quantify actual surface use intensity, and to capture accompanying fine-scale differences in several environmental parameters that often covary 1 or interact with each other across the urban mosaic. Consequently, under the second framework of urbanisation, we used the percentage of ISA surrounding the respective location of every nest analysed in the study. We investigated the impact of urbanisation on alpha diversity as measured by Shannon's diversity, the number of observed OTUs, and Faith's phylogenetic diversity index with linear models (LMs). When analysed based on administrative boundaries (urban vs. rural), Faith's phylogenetic diversity index was lower in urban sites (β = − 0.72 ± 0.36, 95% CI [− 1.44 to − 0.01], p = 0.048, Fig. 4a), while the other alpha diversity metrics did not differ between hosts from urban and rural sites (Supplementary Table 1 5b) dissimilarities of the gut microbiota revealed urban-driven shifts in microbial composition. We statistically tested whether microbial communities differ between nestlings originating from urban and rural areas, different study sites, or with changing ISA levels (fitted as a continuous variable) by PERMANOVA models based on Bray-Curtis and the weighted Uni-Frac dissimilarities. According to the model based on Bray-Curtis dissimilarity, microbial composition differed between urban and rural areas (pseudo-F = 1.66, p = 0.044) and sampling sites (pseudo-F = 1.28, p = 0.019), but not with changing ISA percentages (pseudo-F = 1.08, p = 0.284). PERMANOVA based on Weighted UniFrac distances detected differences in microbial composition between urban and rural sites (pseudo-F = 2.58, p = 0.030), but not between sampling sites (pseudo-F = 1.17, p = 0.202) or with changing ISA percentages (pseudo-F = 1.27, p = 0.218). We also tested the homogeneity of group dispersion between urban and rural sites. PERMDISP based on Bray-Curtis dissimilarity showed homogeneous dispersions for urban and rural samples, indicating that the significant PERMANOVA results were not caused by differences in dispersion among the groups. However, PERMDISP based on weighted UniFrac distances revealed that microbial variation in rural sites is significantly higher than in urban sites (the average distance from the centroids was 0.144 ± 0.006 and 0.125 ± 0.004 and for rural and urban areas, respectively, p = 0.035). Therefore, the significant between-group differences detected by PERMANOVA might reflect differences in location, dispersion, or a combination of the two. There was no significant impact of sampling time in our PERMANOVA model (based on Bray-Curtis dissimilarity: p = 0.30; based on Weighted UniFrac: p = 0.12).
The stacked bar plot revealed prominent taxonomic and compositional differences between urban and rural hosts (Fig. 6a). To determine differentially abundant OTUs between the urban and rural territories, we employed DESeq2 analysis 22 . Overall, we found 20 differentially abundant OTUs (Fig. 6b, Supplementary Table 4). Of these, ten were significantly more abundant in the rural territories, and ten were significantly more abundant in urban territories. Urban hosts exhibited higher abundances of an OTU belonging to the potentially pathogenic microbial family, Enterobacteriaceae 23,24 . Interactive effects of spatial and environmental parameters on microbial change. While nesting sites are here referred to as located in urban or rural sites or areas with varying degrees of ISA, each nestbox is also embedded in a complex web of specific environmental and spatial parameters characterising the urban www.nature.com/scientificreports/ mosaic; these parameters are defined herein as the distance to city centre, closest road and closest path, light pollution, sound pollution, human presence, temperature, tree cover density, and NDVI. These environmental parameters differed between the urban and rural sites (except for the distance to the closest path, NDVI and tree cover density; Supplementary Table 5), and were correlated with ISA (Supplementary Table 3). We investigated which environmental and spatial variables best predicted the microbial alpha diversity (diversity within samples) as measured by Shannon's diversity, the number of OTUs, and Faith's phylogenetic diversity index. As multiple environmental and spatial variables were correlated with each other (Fig. 3b, Supplementary  Table 3), we first examined multicollinearity among all variables and sequentially excluded predictors based on VIF values and biological relevance when analysing predictors of alpha diversity. Our final models included five predictors: distance to the city centre, human presence, temperature, tree cover density and distance to the closest path (Supplementary Table 6). Among these predictors, only tree cover density was positively associated with Shannon's diversity (LMM, β = 0.18 ± 0.09, 95% CI [0.00-0.36], p = 0.047), the number of observed OTUs (β = 1.2 ± 0.59, 95% CI [0.03 to 2.37], p = 0.045) and Faith's phylogenetic diversity (β = 0.56 ± 0.28, 95% CI [0.01 to 1.11], p = 0.048) (Fig. 7).
Before identifying the strongest predictors of between-group microbial dissimilarities, we assessed the spatial autocorrelation of beta diversity by Mantel test. The distance matrix of the geographical coordinates of the sampling sites was correlated with Bray-Curtis dissimilarity matrix (R = 0.054, p = 0.041), but not with the weighted www.nature.com/scientificreports/ UniFrac dissimilarity matrix (R = 0.018, p = 0.250). We further tested the associations between all environmental variables and beta diversity by employing different analytical methods to obtain a coherent vision of the urbandriven variation in gut microbiota. We applied a partial mantel test to investigate the correlations between microbial dissimilarity matrices and the distance matrices of the environmental variables, controlling for the effect of the distance matrix of the geographical coordinates of the sampling sites. We found that sound pollution, distance  www.nature.com/scientificreports/ to the city center, distance to the closest road and tree cover density were significantly correlated with Bray-Curtis and Weighted UniFrac dissimilarities (Table 1). All other parameters were found to be nonsignificant in the partial Mantel tests (Table 1). Based on BIOENV, the best model with the strongest relationship with Bray-Curtis dissimilarity matrix (R = 0.209) similarly contained sound pollution, the distance to the city centre, the distance to the closest road and light pollution (Table 1). However, the best BIOENV model based on weighted UniFrac distance (R = 0.200) contained only sound pollution and distance to the city centre. According to Envfit, patterns observed in the Bray-Curtis and Weighted UniFrac ordination plot were significantly associated with sound pollution, the distance to the city centre, and distance to the closest road, light pollution, tree cover density, and the distance to the closest path,) but not with temperature, NDVI or human presence (Table 1). Overall, there was coherent and conclusive support across all models for sound pollution and distance to the city centre acting as the most important drivers of juvenile gut microbiota variation across the urban mosaic (Fig. 5).

Discussion
Cavity type influences microbial diversity but not community composition. Our dataset, albeit small, indicated that Shannon's diversity, one of the alpha diversity metrics investigated in the study, was higher in nestlings reared in artificial nestboxes, relative to those reared in natural cavities. Since the nestboxes and natural cavities located on a homogenous urban forest (see 18 ), the observed changes cannot be justified by differences in food sources or nest material. Human-made nestboxes are frequently used in ecological studies, as they are convenient substitutes for natural cavities. However, they differ from natural cavities in various aspects, exhibiting very different humidity ranges relative to natural cavities and providing poorer insulation, leading to prominent fluctuations in daily ambient temperatures across 24 h 19 . While the gastrointestinal ecosystem is  www.nature.com/scientificreports/ relatively well-regulated and is not expected to exhibit prominent alterations in response to moderate fluctuations in ambient conditions in endothermic animals, young nestlings of altricial passerines cannot regulate their body temperatures 25 . Therefore, even moderate fluctuations in microclimatic conditions can trigger physiologic stress in nestlings and consequently, influence the immune system 26 . Although we still know very little about how these alterations translate into gut microbiota 27 , some studies revealed an increased microbial diversity with elevated temperatures 28 , probably due to suppression of immune response. Another explanation is that such fluctuations can influence nestling development and differential growth rates between nestlings from different cavity types can ultimately affect the microbial assembly. At the same time, there were no difference in the species richness, phylogenetic diversity, or composition of microbial communities of the nestlings reared in nestboxes and natural cavities. It is important to note that the novel finding regarding increased microbial diversity in nestboxes relies on a small dataset and should be confirmed by larger datasets.
Alterations in the diversity and taxonomic composition of gut microbiota in the urban space. Gut microbial communities differed across the urban mosaic, and we demonstrate that the type of framework used to define the urban space impacts analytical outcomes and results interpretation. Remarkably, occupying urban or rural habitats per se (strictly based on urban/rural administrative definitions) does little to explain changes in alpha diversity. Instead, microbial diversity was found to be lower in highly urbanised space, and was correlated with actual land-use intensity, as measured by ISA (Fig. 4). Thus, the observed alterations in alpha diversity are probably more strongly associated with fine-scale variation in the local habitat, which shows remarkable heterogeneity within the urban mosaic, rather than with factors broadly related to human activities or the geographical locations of the nests. Thus, urban spatial heterogeneity can influence the gut microbiota diversity via different mechanisms (see below). Earlier studies investigating how microbial alpha diversity varies in the urban space have revealed contradictory results. For example, the gut microbial diversity was higher in urban populations of white-crowned sparrows than in their rural counterparts, yet the change was not strongly associated with impervious cover 14,15 , which is inconsistent with our findings. In contrast, adult house sparrows had reduced alpha diversity in urbanised habitats with more than 10% built-up areas 11,29 . Similarly, herring gulls exhibited reduced diversity in highly urbanised areas as determined by the human population 13 . Furthermore, alpha diversity did not differ between urban and rural populations in juvenile house sparrows 29 or was not associated with the percentage of urban land cover in American white ibises 12 . Based on the results reported in this study, we posit that these discrepancies may be driven by differences in the criteria used to define the urban space. Alternatively, but not exclusively, biological differences originating from the taxonomy and life-history stage of the host may also play a role in these discrepancies. For example, avian species with higher tolerance to anthropogenic food might exhibit higher alpha diversity in urban areas. We demonstrate that microbial community composition changes with contrasting land-use patterns. Importantly, and in contrast to the results obtained for alpha diversity, these alterations were prominent only under the urban/rural distinction reflected by administrative borders, indicating that differences in community composition were largely determined by whether the nestling lived in a densely populated urban site. Strikingly, the urban hosts were enriched by a potentially pathogenic microbial family, Enterobacteriaceae 23,30 . A higher prevalence of Enterobacteriaceae has been associated with dysbiosis in mice 31 and higher mortality rates in ostriches 32 . Similarly, urban populations of American white ibises were found to have lower microbial diversity and were more susceptible to Salmonella infections than their rural counterparts 12 . Collectively, these findings indicate that pathogen susceptibility increases in densely populated urban sites, potentially adversely influencing the overall health and fitness of animal hosts. However, it is important to note that the differential abundance analysis did not account for the nonindependence of the samples collected from the same site. Therefore, our results might be related to the increased prevalence of this pathogenic taxon in particular sites rather than overall urban populations.
Taken together, our findings revealed that diversity and composition of the microbial communities of juvenile great tits change across an urban mosaic. Furthermore, using different frameworks to describe the urban space leads to the capturing of different aspects of variation in microbial communities. The initial microbial colonies inhabited by animal hosts in early life are not only critical for the establishment of healthy gut microbiota but also have crucial functions, such as the programming of the immune system 33 , and are involved in the maturation of the nervous system 34 . Therefore, the observed alterations might have long-term consequences on the survival and fitness of hosts 16 . Environmental factors associated with changes in the urban gut microbiota. One of the primary goals of this study was to determine whether distinct environmental variables characterising the urban space can predict changes in the gut microbiota of great tit nestlings. Based on our models, tree cover density was the only environmental axis that covaried with alpha diversity (whilst also confirming that alpha diversity was not spatially or temporally autocorrelated). In line with our findings, white-crowned sparrows occupying territories with greater tree cover also exhibited more diverse gut microbial communities than those occupying regions with less tree cover 14,15 . This is likely to be explained by the tree cover density influencing the diversity and abundance of invertebrates 35 , which are the primary food source for great tits during the chick-rearing period 36 . Reductions in tree cover density, often observed in urban spaces, can also direct animals to search for anthropogenic food 37 . Consistently, the camera recordings of the nests sampled in this study revealed that the amount of anthropogenic food brought into the nests was higher in territories characterized by high ISA percentages than in low-ISA territories (Corsini et al., in prep.), and dietary alterations are known to be associated with changes in alpha diversity 29,38,39 . However, our predictions regarding their diets remain speculative, as we did not analyse the diets of the studied nestlings. Furthermore, it is important to note that tree cover density is correlated with all other variables used in the study (Fig. 3b, Supplementary Table 3 www.nature.com/scientificreports/ to the closest road, NDVI, light and sound pollution had relatively strong associations with tree cover density. Although these variables were not retained in the final models due to multicollinearity issues, the observed alpha diversity patterns might stem from underlying relationships between these variables. An important aim of the study was to determine the environmental parameters predicting compositional changes in gut microbiota. Although seasonal fluctuations can influence gut microbiota 40 , we did not find evidence that temporal variation in sampling leads to compositional changes. However, microbial community composition exhibited spatial structuring. One potential explanation is that the microbial communities residing in the guts of studied nestlings, at least to some extent, represent the environmental pool of microorganisms 41 , which exhibits microgeographic variations 42 . Alternatively, these changes might be related to spatial variation in environmental variables. Based on our models, the two variables that showed the strongest associations with change in the microbial community composition were sound pollution and the distance to the city centre. Sound pollution can influence gut microbiota by different but non-exclusive mechanisms. Noise might impair parent-offspring communication by masking the begging calls of the young, and parents might adjust their food provisioning accordingly 43 . Therefore, the amount of food received by the young might vary depending on the background noise levels, which in turn could influence gut microbiota. Alternatively, acoustic stress can activate the hypothalamic-pituitary-adrenal axis and alter the endocrine profiles 9,44 , causing cascading impacts on glucose metabolism 45 , immune function 46 , oxidative stress 47 and gut physiology 48 and consequently influencing gut microbiota 45,48 . The second strongest predictive variable was the distance to the city center. This interaction might be mediated by environmental variables distributed linearly from the city centre such as light pollution. Indeed, light pollution was a significant predictor based on most of our models and can indirectly affect gut microbiota by altering the host physiology 49 . Alternatively, the observed correlations between the distance to the city center and community composition might be associated through unmeasured chemical pollutants showing collinearity with distance from the city centre, such as particulate matter in ambient air or trace metals 50 . Indeed, environmental chemical exposures have been shown to elicit negative effects on avian physiology 51 and consequently can influence gut microbiota 52 . It is important to note that due to the multicollinearity among these variables, it is not possible to infer the exact correlative factors explaining changes in the microbial community composition, and each of these assumptions requires further experimental testing to infer the exact causative factors driving the observed differences in beta diversity.

Conclusions and outlook
We quantified the extent to which anthropogenic modifications pertaining to the urban environment -including the use of artificial nestboxes-reshapes animal-microbe interactions in great tit nestlings. We report clear trends regarding urban-related shifts in taxonomic composition, a reduction in alpha diversity, and increased pathogen susceptibility. While the exact causal interactions between the environmental parameters and observed microbial changes remain to be tested experimentally, we outlined different mechanisms that are likely to be the driving forces of the measured microbiota variation in the studied urban space. This study thus fills an important gap in providing pioneering evidence on how multi-dimensional anthropogenic changes may relate to gut microbiota during early nestling development, the most critical life-history stage during which gut microbial communities are established and program several developmental processes of their hosts 16,34,53 . Another novel piece of information provided by our study is how the use of nestboxes, the golden standard of ecological field studies, can affect the outcome of research investigating gut microbial communities, an important point that should be considered in future investigations. Our study further adds to the existing evidence that how the urban space is defined might have a prominent impact on the outcome of ecological studies.
Although our study shed light on how anthropogenic change can influence gut microbiota assembly in early life, it also has some potential limitations. First, our relatively small sample size did not allow us to use more sophisticated analysis methods such as machine learning implications to gain deeper insights into drivers of microbial change in urban space. Another drawback regarding the sample size is that our conclusions on cavity type's influence on gut microbiota rely on a small dataset. Consequently, these preliminary findings should be interpreted with caution until they are confirmed by larger datasets. Second, our sampling design was somehow unbalanced, with a larger number of samples collected from urban sites. Although unlikely, we cannot rule out the possible effect of unequal sample size on study outcomes. Third, our sampling was not replicated spatially and temporally; ultimately, the results presented in this work should be assessed by future studies conducted in multiple cities and across multiple years.
Importantly, our study provides important evidence that highlights future opportunities and further refinement in research questions. First, it is pivotal to conduct experimental testing in which the confounding variables are strictly controlled to infer the exact causal effects of the identified environmental factors covarying with microbial changes. Second, whether and how gut microbes can facilitate the adaptation of their hosts to changing environments is an important point that remains to be investigated to achieve full comprehension of the longterm fitness consequences of host-microbe interactions. Third, our study raises the possibility of answering a central question in microbial ecology by leveraging the urban-rural gradient; cross-fostering experiments among contrasting habitats would allow us to disentangle the relative importance of host-specific and environmental factors in shaping microbial communities. Answering these questions would allow us to fully comprehend the ecological and evolutionary dynamics and consequences of animal-microbe symbiosis, in the context of globally increasing urbanisation, specifically in the fields of urban planning, wildlife and public health management.

Materials and methods
Ethics declarations. The experimental protocols, including handling the birds, were approved by the Regional Directorate for Environmental Protection (RDOŚ) in Warsaw, Poland. All experiments were carried out in accordance with relevant guidelines and regulations.

Study sites.
Great tits are hole-nesting passerine birds known to occupy a wide range of habitats, from primary forests 54 to urban city centres 55 . In this species, the nesting environment is an important driver that shapes gut microbiota during early development 56 , making the great tit an ideal subject to study microbial changes across an urban mosaic. Here, we monitored the great tit reproductive cycle in the capital city of Warsaw in Poland, where 565 woodcrete Schwegler nestboxes (type 1b) were set up at nine sites located within and outside the city ( Fig. 1; Supplementary Table 5). We included the following sampling sites: a natural forest (c. 10 km from the city borders and approximately 20 km from the city centre), a peri-urban village bordering a natural forest, two residential areas, two urban woodlands, a large urban park and an office area (2019, 2020). Furthermore, natural cavities located on a Primeval Forest within the urban administrative border were also monitored (Fig. 4).
Capturing urbanisation. The urban space was defined with two different frameworks: under the first framework, each nest was categorised depending on whether it was located within the administrative limits of the city (urban) or outside of these limits (rural). Based on this dichotomous categorization, the urban sites included an office area, two residential areas, an urban forest, an urban park and two urban woodlands. The rural sites comprised a natural forest and a peri-urban village (Supplementary Table 5). Under the second framework, we aimed to capture the heterogeneity in actual surface use patterns and defined urbanisation using the percentage of impervious surface area (ISA) surrounding the respective location of every nest analysed in the study. A map of ISA in Warsaw (in a range from 0 to 100%) with a 20-m pixel resolution was downloaded from Copernicus Land Monitoring Services. ISA included all built-up areas and soil-sealing surfaces substituting original or semi-natural surfaces. Averaged ISA values at each nest were obtained using a 100-m-radius buffer starting from each nestbox location and using geographic information system tools (GIS). The radius was selected on literature-based estimates of parental foraging: blue tits travel about 53.2 m (22.9 SD) in food-poor (but natural) environments while feeding the nestlings 57 ; importantly, the same study has shown that birds also flied beyond 50 m from the nest location in c.1/3 of their foraging trips 57 , therefore, a 100 m radius around each nestbox constitutes a conservative estimate of food foraging distances covered by parents while feeding the nestlings. ISA percentage was used as a continuous variable in alpha and beta diversity analyses.
Environmental and spatial variables. We collected ten environmental and spatial variables as described previously 1,58,59 (also see Supplementary Text 1 for details). Regarding the environmental variables collected on the ground, the (1) human presence was derived by quantifying all humans and dogs, with repeated 30-s long counts performed within a 15-m radius around each nestbox (as detailed in Corsini et al. 59 ). (2) Sound pollution was obtained after averaging recordings that occurred on the DbC scale using hand-held sound level meters equipped with a microphone, over four days throughout the field season, three times per day. (3) Temperature was obtained from 22 Thermocrones ibuttons DS1921G set from 24/04/18 until 30/06/18with a 1-h sampling frequency, and distributed across the entire gradient of urbanisation, following Szulkin et al. 1 . Spatial variables such as the distances from each nestbox to (4) the closest road and to (5) the closest path (for vehicular and pedestrian use, respectively) and (6) the city centre (i.e., the location of the Palace of Culture and Sciences) were computed and measured in metres using the "Measure line" tool in QGIS, as described by Corsini et al. 59 . Furthermore, (7) a distance matrix of the geographical coordinates of the sampling sites was calculated based on the Haversine distances.
We also extracted environmental variables using digital photography and satellite imagery techniques and computed the values of these variables at the nest level in a 100-m-radius buffer. For (8) light pollution, a map of light pollution in Warsaw with a 10-m pixel resolution was extrapolated from night-time digital photographic images shot on 08/10/2015 by astronauts from the International Space Station 60 . For (9) tree cover density, a map of tree cover density in Warsaw (in a range from 0 to 100%) with a 20-m pixel resolution was downloaded from Copernicus Land Monitoring Services. For a proxy of live green vegetation, (10) the normalized difference vegetation index (NDVI) was estimated using satellite images derived from SENTINEL2 that are available on the Earth Explorer website (https:// earth explo rer. usgs. gov).
Sample collection, DNA extraction and library preparation. Faecal samples were collected from 15-day-old great tit chicks hatched in 89 nests (eight and 81 nests were located in natural cavities and nestboxes, respectively; Supplementary Table 5), between the 19th of May 2018 and the 20th of July 2018. Faeces of one chick per nest (from the largest chick in the brood) was collected, except for four nests where more than one chick per nest was sampled (all originating from woodcrete nestboxes). The samples were collected noninvasively while handling the chicks and directly deposited in 5-mL sterile Eppendorf tubes, each filled with 3 mL of RNAlater (Qiagen) and further stored at − 20 °C.
Details of the microbial analysis are described in Maraci et al. 61 . In short, microbial DNA was extracted using the QIAamp PowerFecal DNA Kit (Qiagen), as described in the manufacturer's protocol. The hypervariable V3-V4 region of the 16S ribosomal RNA (rRNA) gene was targeted following the Illumina 16S Metagenomic Library Preparation Guide, 15044223-B. The final amplicon pool contained a pool of blank controls for DNA extraction and PCR amplification and one replicate of a single sample, alongside 94 biological samples sequenced on the Illumina MiSeq system (Illumina, Inc., San Diego, CA, USA). www.nature.com/scientificreports/ Bioinformatic processing. The bioinformatic processing was carried out as described in detail 62 . In short, MiSeq PE reads were assembled in an iterative manner using Flash v1.2.11 63 . All other bioinformatic steps consisted of (1) adapter clipping with cutadapt v1.18 64  Statistical analyses. All consecutive statistical analyses were conducted in R version 4.0.0 68 and Primer-e software version 7 69 . As an initial quality filtering step, all OTUs that could not be classified at the phylum level or that were classified as mitochondria or chloroplasts were discarded. Next, all OTUs that were not represented by at least one read in 2% of the samples or with fewer read counts than 0.001% of the total number of reads were excluded. Subsequently, samples with a lower read count than 900 (N = 12) were removed from the dataset.
We evaluated whether our sample size is sufficient for an accurate estimation using OTU accumulation curve.
To test whether the microbial communities were affected by cavity type, we created a subset of data including samples collected from nestboxes (N = 7) and those collected from natural cavities (N = 6) located in the same urban forest, a highly homogenous ecological site (see 18 , for a map for the distribution of the natural cavities and nestboxes across the study site). Subsequently, the six samples obtained from natural cavities were excluded from the remaining analyses, and the dataset used in further analyses ultimately consisted of 76 samples collected from nestboxes (Supplementary Table 5). Using a two-step approach, we analysed the interactions among microbial community differences, urbanisation and environmental change. First, the urban space was defined using two different frameworks based on (1) the often-used administrative border delineation (urban/rural site) and (2) the percentage of impervious surface area (ISA) surrounding each nest, irrespective of the location of the nest within or outside of the city borders. Thus, microbial alpha and beta diversity was compared (1) between urban and rural sites, and also (2) examined in relation to ISA percentage.
Before the alpha diversity analyses, we rarefied OTU read count data to the lowest read count observed in the dataset (937). Then, we calculated Shannon's diversity index, which accounts for both the abundance and evenness of the taxa present 70 , the number of observed OTUs, which estimates species richness, and Faith's phylogenetic diversity, which incorporates phylogenetic relationships between microbial taxa 71 , as the microbial alpha diversity estimates. We normalised these indexes employing square root transformation. We examined whether these alpha diversity indexes differed between nestboxes and natural cavities by fitting linear models, implemented using the lm function of R package stats 68 .To analyse the impact of urbanisation on alpha diversity, we employed different linear models (LMs) using square-root-transformed Shannon's diversity, the number of observed OTUs, or Faith's phylogenetic diversity indices of the samples collected from nestboxes as the response variable. Depending on the framework used to define urbanisation, we used the urban/rural categorisation or ISA (as a continuous variable) as the fixed effect. Subsequently, to account for the potential impact of differences in the sample collection dates on alpha diversity, we fitted LMs using the number of days between initiation of the study and sample collection as a fixed effect, and Shannon's diversity index, the number of observed OTUs or Faith's phylogenetic diversity as the response variable. The associations between distinct environmental and spatial variables and microbial alpha diversity were also inferred. As an initial step, the extent to which the environmental parameters varied between urban and rural sites was tested with Welch's two-sample t-test. Correlations between alpha diversity, ISA and other environmental variables were then investigated using Pearson's correlation. To analyse the interactions between alpha diversity and environmental variables, three linear mixed models (LMMs) were fitted using Shannon's diversity, the number of observed OTUs, or Faith's phylogenetic diversity indices as the response variables and all environmental variables as fixed effects, as implemented using lme4 version 1.1-25 72 . To examine multicollinearity between the variables, the variance inflation factor (VIF) was calculated with the Car package 73 , and the predictors were sequentially removed from the model based on VIF values and biological relevance, and VIF values were recalculated 74 . This procedure was repeated until all the VIF values were smaller than two. The final models were fitted with the remaining environmental variables as fixed effects, that is: distance to the city centre, human presence, temperature, tree cover density and distance to the closest path. The sampling site was included in these LMMs as random effects to account for the nonindependence of the samples originating from the same sampling site. Furthermore, we also tested for spatial autocorrelation in the simulated scaled residuals of the fitted LMs and LLMs by Moran's I test as implemented in the DHARMa package 75 .
The taxonomic and compositional structures of the microbial communities collected from different localities were visualized with stacked bar plots based on family-level taxonomy using ggplot2 version 3.3.2 76 . To identify differentially abundant OTUs in samples collected from urban and rural localities, the logarithmic fold changes between groups were estimated using a negative binomial Wald test implemented in Deseq2 version 1.12.4 extension 22 of the Phyloseq package version 1.32.0 77 . The significance threshold of the p values was set as 0.05 after a Benjamini and Hochberg false-discovery rate correction 78 .
To analyse community composition 79 , first, the filtered absolute abundance data were normalized by applying cumulative sum scaling (CSS) normalization in the R package Metagenomeseq version 1.30.0 80 to account for unequal sequence coverage. The differences in community composition between urban and rural sites were visualized using non-metric multidimensional scaling (nMDS) based on Bray-Curtis 81 , and the weighted UniFrac 82 dissimilarities implemented using the Vegan package version 2.5 83 . We examined microbial composition in relation to whether samples originated from urban or rural areas, sampling site, ISA percentage (as a continuous value) and sampling time using permutational multivariate analyses of variance (PERMANOVA 84 , with Bray-Curtis and the weighted UniFrac dissimilarities, implemented in Primer-e software with 9999 permutations. Subsequently, using the same dissimilarity measures, we examined the effect of cavity type on microbiota composition by PERMANOVA. We also tested homogeneity of group dispersions using PERMDISP in Primer-e. www.nature.com/scientificreports/ Finally, the interactions between spatial and environmental variables and beta diversity were also evaluated by employing different methods. First, to assess the spatial autocorrelation of beta diversity, we tested the correlations between the Bray-Curtis and weighted UniFrac dissimilarity matrices and the distance matrix of the geographical coordinates of the sampling sites, using the mantel test. Second, we employed a partial mantel test to investigate the correlations between Bray-Curtis and weighted UniFrac microbial dissimilarity matrices and the distance matrices of the environmental variables obtained based on Euclidean distances, controlling for the effect of the distance matrix of the geographical coordinates of the sampling sites. Third, to identify the variables that were strongly related to the first two ordination axes, the multiple regression between each environmental variable and the ordination axes was evaluated using the Envfit function, and the significance of each correlation was tested based on 9999 permutations. Finally, the subset of environmental variables whose Euclidean distance matrices correlated maximally with the microbial distance matrix based on Bray-Curtis was identified employing the BIOENV procedure 85 .

Data availability
Sequences generated for this study have been uploaded to the European Nucleotide Archive (ENA) repository under the accession number: PRJEB44290. The scripts used for processing the data are provided in the GitHub repository at: https:// github. com/ AnnaA ntona touPap/ Urban-relat ed-chang es-in-tit-micro biota. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.